library(tidyverse)
library(readstata13)
library(foreign)
library(readxl)
library(writexl)
library(lubridate)
library(PanelMatch)

# Set working directory
setwd("/Volumes/Backup Plus/Storage-OrganizedFiles/Article1-ReplicationMaterials_FullReplication/1.MainAnalysis")

# Charge the dataset
data1 <- read.dta13("Ortega(2024)-DataAnalysis_Updated(09302024).dta")

data_final <- as.data.frame(data1)
data_final$year <- as.integer(data_final$year)
data_final$divipola <- as.integer(data_final$divipola)

#################################################################################
# Analysis
#################################################################################

# IV: Protests
# A: Create the Matched Set (M)

pm.none5.p <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                         treatment = "iv2.p_strat", refinement.method = "none",
                         data = data_final, match.missing = TRUE,
                         covs.formula = ~ I(lag(cv1.lpop, 1:3)) + I(lag(cv2.rur_p, 1:3)) 
                         + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                         + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                         + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.2.no_p_str, 1:3)) 
                         + I(lag(dv4b.2.1.farc_ckill_up1_r, 1:3)),
                         size.match = 5, qoi = "att", outcome.var = "dv4b.2.1.farc_ckill_up1_r",
                         lead = 0:4, forbid.treatment.reversal = FALSE,
                         use.diagonal.variance.matrix = TRUE)

pm.maha5b.p <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                          treatment = "iv2.p_strat", refinement.method = "mahalanobis",
                          data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                          covs.formula = ~ I(lag(cv1.lpop, 1:3)) + I(lag(cv2.rur_p, 1:3)) 
                          + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                          + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                          + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.2.no_p_str, 1:3)) 
                          + I(lag(dv4b.2.1.farc_ckill_up1_r, 1:3)),
                          size.match = 5, qoi = "att" , outcome.var = "dv4b.2.1.farc_ckill_up1_r",
                          lead = 0:4, forbid.treatment.reversal = FALSE,
                          use.diagonal.variance.matrix = TRUE)

pm.ps.m5.p <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                         treatment = "iv2.p_strat", refinement.method = "ps.match",
                         data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                         covs.formula = ~ I(lag(cv1.lpop, 1:3)) 
                         + I(lag(cv2.rur_p, 1:3)) 
                         + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                         + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                         + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.2.no_p_str, 1:3)) 
                         + I(lag(dv4b.2.1.farc_ckill_up1_r, 1:3)),
                         size.match = 5, qoi = "att", outcome.var = "dv4b.2.1.farc_ckill_up1_r",
                         lead = 0:4, forbid.treatment.reversal = FALSE)

pm.ps.w5.p <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                         treatment = "iv2.p_strat", refinement.method = "ps.weight",
                         data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                         covs.formula = ~ I(lag(cv1.lpop, 1:3)) 
                         + I(lag(cv2.rur_p, 1:3)) 
                         + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                         + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                         + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.2.no_p_str, 1:3)) 
                         + I(lag(dv4b.2.1.farc_ckill_up1_r, 1:3)),
                         size.match = 5, qoi = "att", outcome.var = "dv4b.2.1.farc_ckill_up1_r",
                         lead = 0:4, forbid.treatment.reversal = FALSE)

pm.CBPSm5.p <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                          treatment = "iv2.p_strat", refinement.method = "CBPS.match",
                          data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                          covs.formula = ~ I(lag(cv1.lpop, 1:3)) + I(lag(cv2.rur_p, 1:3)) 
                          + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                          + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                          + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.2.no_p_str, 1:3)) 
                          + I(lag(dv4b.2.1.farc_ckill_up1_r, 1:3)),
                          size.match = 5, qoi = "att" , outcome.var = "dv4b.2.1.farc_ckill_up1_r",
                          lead = 0:4, forbid.treatment.reversal = FALSE,
                          use.diagonal.variance.matrix = TRUE)

pm.CBPSw5.p <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                          treatment = "iv2.p_strat", refinement.method = "CBPS.weight",
                          data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                          covs.formula = ~ I(lag(cv1.lpop, 1:3)) 
                          + I(lag(cv2.rur_p, 1:3)) 
                          + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                          + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                          + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.2.no_p_str, 1:3)) 
                          + I(lag(dv4b.2.1.farc_ckill_up1_r, 1:3)),
                          size.match = 5, qoi = "att", outcome.var = "dv4b.2.1.farc_ckill_up1_r",
                          lead = 0:4, forbid.treatment.reversal = FALSE)

# IV: Violent Self-Protection
# A: Create the Matched Set (M)
pm.none5.v <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                         treatment = "iv4b.v_strat", refinement.method = "none",
                         data = data_final, match.missing = TRUE,
                         covs.formula = ~ I(lag(cv1.lpop, 1:3)) + I(lag(cv2.rur_p, 1:3)) 
                         + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                         + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                         + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.4b.no_v_str, 1:3)) 
                         + I(lag(dv4b.2.1.farc_ckill_up1_r, 1:3)),
                         size.match = 5, qoi = "att", outcome.var = "dv4b.2.1.farc_ckill_up1_r",
                         lead = 0:4, forbid.treatment.reversal = FALSE,
                         use.diagonal.variance.matrix = TRUE)

pm.maha5b.v <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                          treatment = "iv4b.v_strat", refinement.method = "mahalanobis",
                          data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                          covs.formula = ~ I(lag(cv1.lpop, 1:3)) + I(lag(cv2.rur_p, 1:3)) 
                          + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                          + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                          + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.4b.no_v_str, 1:3)) 
                          + I(lag(dv4b.2.1.farc_ckill_up1_r, 1:3)),
                          size.match = 5, qoi = "att" , outcome.var = "dv4b.2.1.farc_ckill_up1_r",
                          lead = 0:4, forbid.treatment.reversal = FALSE,
                          use.diagonal.variance.matrix = TRUE)

pm.ps.m5.v <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                         treatment = "iv4b.v_strat", refinement.method = "ps.match",
                         data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                         covs.formula = ~ I(lag(cv1.lpop, 1:3)) 
                         + I(lag(cv2.rur_p, 1:3)) 
                         + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                         + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                         + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.4b.no_v_str, 1:3)) 
                         + I(lag(dv4b.2.1.farc_ckill_up1_r, 1:3)),
                         size.match = 5, qoi = "att", outcome.var = "dv4b.2.1.farc_ckill_up1_r",
                         lead = 0:4, forbid.treatment.reversal = FALSE)

pm.ps.w5.v <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                         treatment = "iv4b.v_strat", refinement.method = "ps.weight",
                         data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                         covs.formula = ~ I(lag(cv1.lpop, 1:3)) 
                         + I(lag(cv2.rur_p, 1:3)) 
                         + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                         + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                         + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.4b.no_v_str, 1:3)) 
                         + I(lag(dv4b.2.1.farc_ckill_up1_r, 1:3)),
                         size.match = 5, qoi = "att", outcome.var = "dv4b.2.1.farc_ckill_up1_r",
                         lead = 0:4, forbid.treatment.reversal = FALSE)

pm.CBPSm5.v <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                          treatment = "iv4b.v_strat", refinement.method = "CBPS.match",
                          data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                          covs.formula = ~ I(lag(cv1.lpop, 1:3)) + I(lag(cv2.rur_p, 1:3)) 
                          + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                          + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                          + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.4b.no_v_str, 1:3)) 
                          + I(lag(dv4b.2.1.farc_ckill_up1_r, 1:3)),
                          size.match = 5, qoi = "att" , outcome.var = "dv4b.2.1.farc_ckill_up1_r",
                          lead = 0:4, forbid.treatment.reversal = FALSE,
                          use.diagonal.variance.matrix = TRUE)

pm.CBPSw5.v <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                          treatment = "iv4b.v_strat", refinement.method = "CBPS.weight",
                          data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                          covs.formula = ~ I(lag(cv1.lpop, 1:3)) 
                          + I(lag(cv2.rur_p, 1:3)) 
                          + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                          + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                          + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.4b.no_v_str, 1:3)) 
                          + I(lag(dv4b.2.1.farc_ckill_up1_r, 1:3)),
                          size.match = 5, qoi = "att", outcome.var = "dv4b.2.1.farc_ckill_up1_r",
                          lead = 0:4, forbid.treatment.reversal = FALSE)

# IV: Sanctuary
# A: Create the Matched Set (M)
pm.none5.s <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                         treatment = "iv3.s_strat", refinement.method = "none",
                         data = data_final, match.missing = TRUE,
                         covs.formula = ~ I(lag(cv1.lpop, 1:3)) + I(lag(cv2.rur_p, 1:3)) 
                         + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                         + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                         + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.3.no_s_str, 1:3)) 
                         + I(lag(dv4b.2.1.farc_ckill_up1_r, 1:3)),
                         size.match = 5, qoi = "att", outcome.var = "dv4b.2.1.farc_ckill_up1_r",
                         lead = 0:4, forbid.treatment.reversal = FALSE,
                         use.diagonal.variance.matrix = TRUE)

pm.maha5b.s <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                          treatment = "iv3.s_strat", refinement.method = "mahalanobis",
                          data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                          covs.formula = ~ I(lag(cv1.lpop, 1:3)) + I(lag(cv2.rur_p, 1:3)) 
                          + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                          + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                          + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.3.no_s_str, 1:3)) 
                          + I(lag(dv4b.2.1.farc_ckill_up1_r, 1:3)),
                          size.match = 5, qoi = "att" , outcome.var = "dv4b.2.1.farc_ckill_up1_r",
                          lead = 0:4, forbid.treatment.reversal = FALSE,
                          use.diagonal.variance.matrix = TRUE)

pm.ps.m5.s <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                         treatment = "iv3.s_strat", refinement.method = "ps.match",
                         data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                         covs.formula = ~ I(lag(cv1.lpop, 1:3)) 
                         + I(lag(cv2.rur_p, 1:3)) 
                         + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                         + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                         + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.3.no_s_str, 1:3)) 
                         + I(lag(dv4b.2.1.farc_ckill_up1_r, 1:3)),
                         size.match = 5, qoi = "att", outcome.var = "dv4b.2.1.farc_ckill_up1_r",
                         lead = 0:4, forbid.treatment.reversal = FALSE)

pm.ps.w5.s <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                         treatment = "iv3.s_strat", refinement.method = "ps.weight",
                         data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                         covs.formula = ~ I(lag(cv1.lpop, 1:3)) 
                         + I(lag(cv2.rur_p, 1:3)) 
                         + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                         + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                         + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.3.no_s_str, 1:3)) 
                         + I(lag(dv4b.2.1.farc_ckill_up1_r, 1:3)),
                         size.match = 5, qoi = "att", outcome.var = "dv4b.2.1.farc_ckill_up1_r",
                         lead = 0:4, forbid.treatment.reversal = FALSE)

pm.CBPSm5.s <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                          treatment = "iv3.s_strat", refinement.method = "CBPS.match",
                          data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                          covs.formula = ~ I(lag(cv1.lpop, 1:3)) + I(lag(cv2.rur_p, 1:3)) 
                          + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                          + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                          + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.3.no_s_str, 1:3)) 
                          + I(lag(dv4b.2.1.farc_ckill_up1_r, 1:3)),
                          size.match = 5, qoi = "att" , outcome.var = "dv4b.2.1.farc_ckill_up1_r",
                          lead = 0:4, forbid.treatment.reversal = FALSE,
                          use.diagonal.variance.matrix = TRUE)

pm.CBPSw5.s <- PanelMatch(lag = 3, time.id = "year", unit.id = "divipola",
                          treatment = "iv3.s_strat", refinement.method = "CBPS.weight",
                          data = data_final, match.missing = FALSE, listwise.delete = TRUE,
                          covs.formula = ~ I(lag(cv1.lpop, 1:3)) 
                          + I(lag(cv2.rur_p, 1:3)) 
                          + I(lag(cv3c.c_lsh1, 1:3)) + I(lag(cv4.nbi, 1:3))
                          + I(lag(cv5c.1.bd_o, 1:3)) + I(lag(cv6b.1.ce_o_r, 1:3))
                          + I(lag(cv7b.3.pgf_cd_o_r, 1:3)) + I(lag(cv8.3.no_s_str, 1:3)) 
                          + I(lag(dv4b.2.1.farc_ckill_up1_r, 1:3)),
                          size.match = 5, qoi = "att", outcome.var = "dv4b.2.1.farc_ckill_up1_r",
                          lead = 0:4, forbid.treatment.reversal = FALSE)

#### Plots ####

msets.none5.p <- pm.none5.p$att
msets.maha5b.p <- pm.maha5b.p$att
msets.ps.m5.p <- pm.ps.m5.p$att
msets.ps.w5.p <- pm.ps.w5.p$att
msets.CBPSm5.p <- pm.CBPSm5.p$att
msets.CBPSw5.p <- pm.CBPSw5.p$att

# C Panel Estimate
pe.results5.p2_maha5b <- PanelEstimate(sets = pm.maha5b.p, data = data_final,
                                       se.method = "conditional",
                                       number.iterations = 1000,
                                       confidence.level = .95)

pe.results5.p2_ps.m5 <- PanelEstimate(sets = pm.ps.m5.p, data = data_final,
                                      se.method = "conditional",
                                      number.iterations = 1000,
                                      confidence.level = .95)

pe.results5.p2_ps.w5 <- PanelEstimate(sets = pm.ps.w5.p, data = data_final,
                                      se.method = "conditional",
                                      number.iterations = 1000,
                                      confidence.level = .95)

pe.results5.p2_CBPSm5 <- PanelEstimate(sets = pm.CBPSm5.p, data = data_final,
                                       se.method = "conditional",
                                       number.iterations = 1000,
                                       confidence.level = .95)

pe.results5.p2_CBPSw5 <- PanelEstimate(sets = pm.CBPSw5.p, data = data_final,
                                       se.method = "conditional",
                                       number.iterations = 1000,
                                       confidence.level = .95)

msets.none5.v <- pm.none5.v$att
msets.maha5b.v <- pm.maha5b.v$att
msets.ps.m5.v <- pm.ps.m5.v$att
msets.ps.w5.v <- pm.ps.w5.v$att
msets.CBPSm5.v <- pm.CBPSm5.v$att
msets.CBPSw5.v <- pm.CBPSw5.v$att

pe.results5.v2_maha5b <- PanelEstimate(sets = pm.maha5b.v, data = data_final,
                                       se.method = "conditional",
                                       number.iterations = 1000,
                                       confidence.level = .95)

pe.results5.v2_ps.m5 <- PanelEstimate(sets = pm.ps.m5.v, data = data_final,
                                      se.method = "conditional",
                                      number.iterations = 1000,
                                      confidence.level = .95)

pe.results5.v2_ps.w5 <- PanelEstimate(sets = pm.ps.w5.v, data = data_final,
                                      se.method = "conditional",
                                      number.iterations = 1000,
                                      confidence.level = .95)

pe.results5.v2_CBPSm5 <- PanelEstimate(sets = pm.CBPSm5.v, data = data_final,
                                       se.method = "conditional",
                                       number.iterations = 1000,
                                       confidence.level = .95)

pe.results5.v2_CBPSw5 <- PanelEstimate(sets = pm.CBPSw5.v, data = data_final,
                                       se.method = "conditional",
                                       number.iterations = 1000,
                                       confidence.level = .95)

msets.none5.s <- pm.none5.s$att
msets.maha5b.s <- pm.maha5b.s$att
msets.ps.m5.s <- pm.ps.m5.s$att
msets.ps.w5.s <- pm.ps.w5.s$att
msets.CBPSm5.s <- pm.CBPSm5.s$att
msets.CBPSw5.s <- pm.CBPSw5.s$att

pe.results5.s2_maha5b <- PanelEstimate(sets = pm.maha5b.s, data = data_final,
                                       se.method = "conditional",
                                       number.iterations = 1000,
                                       confidence.level = .95)

pe.results5.s2_ps.m5 <- PanelEstimate(sets = pm.ps.m5.s, data = data_final,
                                      se.method = "conditional",
                                      number.iterations = 1000,
                                      confidence.level = .95)

pe.results5.s2_ps.w5 <- PanelEstimate(sets = pm.ps.w5.s, data = data_final,
                                      se.method = "conditional",
                                      number.iterations = 1000,
                                      confidence.level = .95)

pe.results5.s2_CBPSm5 <- PanelEstimate(sets = pm.CBPSm5.s, data = data_final,
                                       se.method = "conditional",
                                       number.iterations = 1000,
                                       confidence.level = .95)

pe.results5.s2_CBPSw5 <- PanelEstimate(sets = pm.CBPSw5.s, data = data_final,
                                       se.method = "conditional",
                                       number.iterations = 1000,
                                       confidence.level = .95)

##################
# Strategic Factors
##################

data_final <- data_final %>%
  mutate(zone4 = case_when(cv5c.1.bd_o < 0.34 ~ "Dominant",
                           TRUE ~ "Nondominant"))

## Protest ##
pe.results5.p2_maha5b_mod1 <- PanelEstimate(sets = pm.maha5b.p, data = data_final, 
                                            moderator = "zone4",
                                            se.method = "conditional",
                                            number.iterations = 1000,
                                            confidence.level = .95)

pe.results5.p2_ps.m5_mod1 <- PanelEstimate(sets = pm.ps.m5.p, data = data_final,
                                           moderator = "zone4",
                                           se.method = "conditional",
                                           number.iterations = 1000,
                                           confidence.level = .95)

pe.results5.p2_ps.w5_mod1 <- PanelEstimate(sets = pm.ps.w5.p, data = data_final,
                                           moderator = "zone4",
                                           se.method = "conditional",
                                           number.iterations = 1000,
                                           confidence.level = .95)

pe.results5.p2_CBPSm5_mod1 <- PanelEstimate(sets = pm.CBPSm5.p, data = data_final,
                                            moderator = "zone4",
                                            se.method = "conditional",
                                            number.iterations = 1000,
                                            confidence.level = .95)

pe.results5.p2_CBPSw5_mod1 <- PanelEstimate(sets = pm.CBPSw5.p, data = data_final,
                                            moderator = "zone4",
                                            se.method = "conditional",
                                            number.iterations = 1000,
                                            confidence.level = .95)

## Sanctuary ##
pe.results5.s2_maha5b_mod1 <- PanelEstimate(sets = pm.maha5b.s, data = data_final,
                                            moderator = "zone4",
                                            se.method = "conditional",
                                            number.iterations = 1000,
                                            confidence.level = .95)

pe.results5.s2_ps.m5_mod1 <- PanelEstimate(sets = pm.ps.m5.s, data = data_final,
                                           moderator = "zone4",
                                           se.method = "conditional",
                                           number.iterations = 1000,
                                           confidence.level = .95)

pe.results5.s2_ps.w5_mod1 <- PanelEstimate(sets = pm.ps.w5.s, data = data_final,
                                           moderator = "zone4",
                                           se.method = "conditional",
                                           number.iterations = 1000,
                                           confidence.level = .95)

pe.results5.s2_CBPSm5_mod1 <- PanelEstimate(sets = pm.CBPSm5.s, data = data_final,
                                            moderator = "zone4",
                                            se.method = "conditional",
                                            number.iterations = 1000,
                                            confidence.level = .95)

pe.results5.s2_CBPSw5_mod1 <- PanelEstimate(sets = pm.CBPSw5.s, data = data_final,
                                            moderator = "zone4",
                                            se.method = "conditional",
                                            number.iterations = 1000,
                                            confidence.level = .95)

## VSP ##

pe.results5.v2_maha5b_mod1 <- PanelEstimate(sets = pm.maha5b.v, data = data_final,
                                            moderator = "zone4",      
                                            se.method = "conditional",
                                            number.iterations = 1000,
                                            confidence.level = .95)

pe.results5.v2_ps.m5_mod1 <- PanelEstimate(sets = pm.ps.m5.v, data = data_final,
                                           moderator = "zone4",      
                                           se.method = "conditional",
                                           number.iterations = 1000,
                                           confidence.level = .95)

pe.results5.v2_ps.w5_mod1 <- PanelEstimate(sets = pm.ps.w5.v, data = data_final,
                                           moderator = "zone4",      
                                           se.method = "conditional",
                                           number.iterations = 1000,
                                           confidence.level = .95)

pe.results5.v2_CBPSm5_mod1 <- PanelEstimate(sets = pm.CBPSm5.v, data = data_final,
                                            moderator = "zone4",      
                                            se.method = "conditional",
                                            number.iterations = 1000,
                                            confidence.level = .95)

pe.results5.v2_CBPSw5_mod1 <- PanelEstimate(sets = pm.CBPSw5.v, data = data_final,
                                            moderator = "zone4",      
                                            se.method = "conditional",
                                            number.iterations = 1000,
                                            confidence.level = .95)

##################
# Organizational Factors I: Ideology
##################

data_final <- data_final %>%
  mutate(l_supp = case_when(cv3c.c_lsh1 > median_c_lsh_m ~ "High",
                            cv3c.c_lsh1 < median_c_lsh_m | cv3c.c_lsh1 == 0 ~ "Low"))

## Protest ##
pe.results5.p2_maha5b_mod2 <- PanelEstimate(sets = pm.maha5b.p, data = data_final, 
                                            moderator = "l_supp",
                                            se.method = "conditional",
                                            number.iterations = 1000,
                                            confidence.level = .95)

pe.results5.p2_ps.m5_mod2 <- PanelEstimate(sets = pm.ps.m5.p, data = data_final,
                                           moderator = "l_supp",
                                           se.method = "conditional",
                                           number.iterations = 1000,
                                           confidence.level = .95)

pe.results5.p2_ps.w5_mod2 <- PanelEstimate(sets = pm.ps.w5.p, data = data_final,
                                           moderator = "l_supp",
                                           se.method = "conditional",
                                           number.iterations = 1000,
                                           confidence.level = .95)

pe.results5.p2_CBPSm5_mod2 <- PanelEstimate(sets = pm.CBPSm5.p, data = data_final,
                                            moderator = "l_supp",
                                            se.method = "conditional",
                                            number.iterations = 1000,
                                            confidence.level = .95)

pe.results5.p2_CBPSw5_mod2 <- PanelEstimate(sets = pm.CBPSw5.p, data = data_final,
                                            moderator = "l_supp",
                                            se.method = "conditional",
                                            number.iterations = 1000,
                                            confidence.level = .95)

## Sanctuary ##
pe.results5.s2_maha5b_mod2 <- PanelEstimate(sets = pm.maha5b.s, data = data_final,
                                            moderator = "l_supp",
                                            se.method = "conditional",
                                            number.iterations = 1000,
                                            confidence.level = .95)

pe.results5.s2_ps.m5_mod2 <- PanelEstimate(sets = pm.ps.m5.s, data = data_final,
                                           moderator = "l_supp",
                                           se.method = "conditional",
                                           number.iterations = 1000,
                                           confidence.level = .95)

pe.results5.s2_ps.w5_mod2 <- PanelEstimate(sets = pm.ps.w5.s, data = data_final,
                                           moderator = "l_supp",
                                           se.method = "conditional",
                                           number.iterations = 1000,
                                           confidence.level = .95)

pe.results5.s2_CBPSm5_mod2 <- PanelEstimate(sets = pm.CBPSm5.s, data = data_final,
                                            moderator = "l_supp",
                                            se.method = "conditional",
                                            number.iterations = 1000,
                                            confidence.level = .95)

pe.results5.s2_CBPSw5_mod2 <- PanelEstimate(sets = pm.CBPSw5.s, data = data_final,
                                            moderator = "l_supp",
                                            se.method = "conditional",
                                            number.iterations = 1000,
                                            confidence.level = .95)

##################
# Organizational Factors II: Recruitment.
##################

data_final <- data_final %>%
  mutate(reb_ca = case_when(reb_ru == 0 ~ "No Recruitment",
                            reb_ru > 0 ~ "Recruitment"))

##
pe.results5.p2_maha5b_mod3 <- PanelEstimate(sets = pm.maha5b.p, data = data_final, 
                                            moderator = "reb_ca",
                                            se.method = "conditional",
                                            number.iterations = 1000,
                                            confidence.level = .95)

pe.results5.p2_ps.m5_mod3 <- PanelEstimate(sets = pm.ps.m5.p, data = data_final,
                                           moderator = "reb_ca",
                                           se.method = "conditional",
                                           number.iterations = 1000,
                                           confidence.level = .95)

pe.results5.p2_ps.w5_mod3 <- PanelEstimate(sets = pm.ps.w5.p, data = data_final,
                                           moderator = "reb_ca",
                                           se.method = "conditional",
                                           number.iterations = 1000,
                                           confidence.level = .95)

pe.results5.p2_CBPSm5_mod3 <- PanelEstimate(sets = pm.CBPSm5.p, data = data_final,
                                            moderator = "reb_ca",
                                            se.method = "conditional",
                                            number.iterations = 1000,
                                            confidence.level = .95)

pe.results5.p2_CBPSw5_mod3 <- PanelEstimate(sets = pm.CBPSw5.p, data = data_final,
                                            moderator = "reb_ca",
                                            se.method = "conditional",
                                            number.iterations = 1000,
                                            confidence.level = .95)

## Sanctuary ##
pe.results5.s2_maha5b_mod3 <- PanelEstimate(sets = pm.maha5b.s, data = data_final,
                                            moderator = "reb_ca",
                                            se.method = "conditional",
                                            number.iterations = 1000,
                                            confidence.level = .95)

pe.results5.s2_ps.m5_mod3 <- PanelEstimate(sets = pm.ps.m5.s, data = data_final,
                                           moderator = "reb_ca",
                                           se.method = "conditional",
                                           number.iterations = 1000,
                                           confidence.level = .95)

pe.results5.s2_ps.w5_mod3 <- PanelEstimate(sets = pm.ps.w5.s, data = data_final,
                                           moderator = "reb_ca",
                                           se.method = "conditional",
                                           number.iterations = 1000,
                                           confidence.level = .95)

pe.results5.s2_CBPSm5_mod3 <- PanelEstimate(sets = pm.CBPSm5.s, data = data_final,
                                            moderator = "reb_ca",
                                            se.method = "conditional",
                                            number.iterations = 1000,
                                            confidence.level = .95)

pe.results5.s2_CBPSw5_mod3 <- PanelEstimate(sets = pm.CBPSw5.s, data = data_final,
                                            moderator = "reb_ca",
                                            se.method = "conditional",
                                            number.iterations = 1000,
                                            confidence.level = .95)


save.image(file='OMC-ResultsMainModel_Updated.RData')

########################################################
## End of this part
########################################################